#library
# library(devtools)
#install_version("mediation", version = "4.0.0", repos = "http://cran.us.r-project.org")
library(mediation)
library(haven)
# library(dummies)
# library(installr)

#dataset
mydata <- read_dta("data/Demo_book_country_3.dta")

# names(mydata)
# baseline model
summary(mod1 <- lm(v2x_polyarchy_imp_100 ~ predallports53_100km  + year, data=mydata))

#mediators: Maddison_gdppc_1990_estimate_ln , nearest_univ_dist_ln_stock_1 , eur_lang_pct , e_micowimp_ln

####GDP

mydata$mediator <- mydata$Maddison_gdppc_1990_estimate_ln

#model of the mediator
med.fit1 <- lm(mediator~ predallports53_100km  + year , data = mydata)

#model of the outcome
out.fit1 <- lm(v2x_polyarchy_imp_100 ~ mediator + predallports53_100km  + year, data = mydata)

med.out2 <- mediate(med.fit1, out.fit1, treat = "predallports53_100km", mediator = "mediator", sims = 1000, dropobs = TRUE)

summary(med.out2)


## tiff("output/mediator1.tiff")
## plot(med.out2, xlim=c(-.004,0), main="Mediator: GDP")
## dev.off()



####UNIVERSITIES

mydata$mediator <- mydata$nearest_univ_dist_ln_stock_1

#model of the mediator
med.fit1 <- lm(mediator~ predallports53_100km  + year , data = mydata)

#model of the outcome
out.fit1 <- lm(v2x_polyarchy_imp_100 ~ mediator + predallports53_100km  + year, data = mydata)

med.out2 <- mediate(med.fit1, out.fit1, treat = "predallports53_100km", mediator = "mediator", sims = 1000, dropobs = TRUE)
summary(med.out2)


## tiff("output/mediator2.tiff")
## plot(med.out2, xlim=c(-.004,0), main="Mediator: Universities")
## dev.off()

####language

mydata$mediator <- mydata$eur_pct_est_smooth_ln

#model of the mediator
med.fit1 <- lm(mediator~ predallports53_100km  + year , data = mydata)

#model of the outcome
out.fit1 <- lm(v2x_polyarchy_imp_100 ~ mediator + predallports53_100km  + year, data = mydata)

med.out2 <- mediate(med.fit1, out.fit1, treat = "predallports53_100km", mediator = "mediator", sims = 1000, dropobs = TRUE)
summary(med.out2)


## tiff("output/mediator3.tiff")
## plot(med.out2, xlim=c(-.004,0), main="Mediator: European language")
## dev.off()

####imports 

mydata$mediator <- mydata$e_micowimp_ln

#model of the mediator
med.fit1 <- lm(mediator~ predallports53_100km  + year , data = mydata)

#model of the outcome
out.fit1 <- lm(v2x_polyarchy_imp_100 ~ mediator + predallports53_100km  + year, data = mydata)

med.out2 <- mediate(med.fit1, out.fit1, treat = "predallports53_100km", mediator = "mediator", sims = 1000, dropobs = TRUE)
summary(med.out2)

## tiff("output/mediator4.tiff")
## plot(med.out2, xlim=c(-.004,0), main="Mediator: Imports")
## dev.off()


####tton milper

mydata$mediator <- mydata$totton_milper

#model of the mediator
med.fit1 <- lm(mediator~ predallports53_100km  + year , data = mydata)

#model of the outcome
out.fit1 <- lm(v2x_polyarchy_imp_100 ~ mediator + predallports53_100km  + year, data = mydata)

med.out2 <- mediate(med.fit1, out.fit1, treat = "predallports53_100km", mediator = "mediator", sims = 1000, dropobs = TRUE)
summary(med.out2)


## tiff("output/mediator5.tiff")
## plot(med.out2, xlim=c(-.004,0), main="Mediator: Naval ship ration")
## dev.off()



#creating a common component factor
fact <- subset(mydata, select=c(e_micowimp_ln, totton_milper, eur_pct_est_smooth_ln , nearest_univ_dist_ln_stock_1, Maddison_gdppc_1990_estimate_ln, predallports53_100km  , year, 
                                v2x_polyarchy_imp_100))
fact <- fact[complete.cases(fact)==T,]
fact2 <- subset(fact, select=c( e_micowimp_ln, totton_milper, eur_pct_est_smooth_ln , nearest_univ_dist_ln_stock_1, Maddison_gdppc_1990_estimate_ln))
fit <- factanal(fact2, 1, rotation="varimax", scores="regression")
fact$mediator <- fit$scores
print(fit)


#model of the mediator
med.fit1 <- lm(mediator~ predallports53_100km  + year , data = fact)

#model of the outcome
out.fit1 <- lm(v2x_polyarchy_imp_100 ~ mediator + predallports53_100km  + year, data = fact)

med.out2 <- mediate(med.fit1, out.fit1, treat = "predallports53_100km", mediator = "mediator", sims = 1000, dropobs = TRUE)
summary(med.out2)

## tiff("output/mediator6.tiff")
## plot(med.out2, xlim=c(-.004,0), main="Mediator: Common factor for all the mediators")
## dev.off()
